Spatio-temporal changes of small protist and free-living bacterial communities in a temperate dimictic lake: insights from metabarcoding and machine learning

Abstract Microbial communities, which include prokaryotes and protists, play an important role in aquatic ecosystems and influence ecological processes. To understand these communities, metabarcoding provides a powerful tool to assess their taxonomic composition and track spatio-temporal dynamics in both marine and freshwater environments. While marine ecosystems have been extensively studied, there is a notable research gap in understanding eukaryotic microbial communities in temperate lakes. Our study addresses this gap by investigating the free-living bacteria and small protist communities in Lake Roś (Poland), a dimictic temperate lake. Metabarcoding analysis revealed that both the bacterial and protist communities exhibit distinct seasonal patterns that are not necessarily shaped by dominant taxa. Furthermore, machine learning and statistical methods identified crucial amplicon sequence variants (ASVs) specific to each season. In addition, we identified a distinct community in the anoxic hypolimnion. We have also shown that the key factors shaping the composition of analysed community are temperature, oxygen, and silicon concentration. Understanding these community structures and the underlying factors is important in the context of climate change potentially impacting mixing patterns and leading to prolonged stratification.


Introduction
Pr otists ar e abundant and div erse eukaryotic micr oor ganisms in aquatic ecosystems and fulfil critical ecosystem functions .T hey play an essential role in organic matter cycling by contributing to primary production and decomposition of organic matter and constitute a link between prokaryotes and higher trophic level or ganisms (Car on 1994, Nakano et al. 1998, Posch et al. 2015, Šimek et al. 2020 ).Despite their ecological significance, the compr ehensiv e understanding of pr otist div ersity r emains limited.Recent advancements in molecular techniques have spurred a surge in diversity studies, revealing an unexpectedly high diversity of protists across various aquatic en vironments , particularly in oceans (e.g.López-García et al. 2001, Lovejoy et al. 2006, Worden et al. 2006, Stoeck et al. 2010, de Vargas et al. 2015, Lima-Mendez et al. 2015, Massana et al. 2015, Seeleuthner et al. 2018, Sunagawa et al. 2020 ).Howe v er, studies on fr eshwater pr otist div ersity r emain compar ativ el y scarce, often focusing on specific types of water bodies (e.g.Charvet et al. 2012, Cruaud et al. 2019, David et al. 2021, Metz et al. 2022 ).Furthermore, small protists in particular, although recognized as important components of micr obial comm unities in lacustrine en vironments (F enchel 1986(F enchel , Stoc kner 1988 ) ), wer e not studied in detail befor e the adv ent of molecular methods .T hese micr oor ganisms ar e often too small to be easily identified and lack distinct morphological features, so their true diversity was inaccessible and their taxonomy poorly understood.
F reshw ater ecosystems are more fragmented and isolated (Dodson 1992, Reche et al. 2005 ), compared to the ocean, where micr obial comm unities ar e disseminated on a global scale via ocean currents (Villarino et al. 2018, Richter et al. 2022 ).This intrinsic lo w er connectivity of freshw ater ecosystems hinders the dispersal of freshwater organisms and increases their genetic diversity (Manel et al. 2020, Miller 2021 ).Furthermor e, fr eshwater ecosystems' environmental conditions are more heterogeneous and m uc h mor e sensitiv e to external factors than those in the oceans (Simon et al. 2015b ).Recent analyses across div erse habitats r e v ealed a ppar ent differ ences in the taxonomic composition of the major protistan lineages and a higher βdiversity in freshwater bodies than in the other systems (Singer et al. 2021, Xiong et al. 2021 ), prompting studies on freshwater ecosystems.
Within the realm of freshwater ecosystems, lakes are the most studied (Charvet et al. 2012, Lepère et al. 2016, Boenigk et al. 2018 ).Notabl y, r esearc h has pr edominantl y concentr ated on high mountain lak es (Filk er et al. 2016, Kammerlander et al. 2016, Boenigk et al. 2018 ) and polar lakes (Daniel et al. 2016, Stoof-Leichsenring et al. 2020 ) due to their extreme conditions, including temper atur e , nutrient a v ailability, and UV r adiation.Se v er al studies have been performed on shallow eutrophic lakes (Simon et al. 2015a ,b ), lak es with ano xic hypolimnion (Oikonomou et al. 2015, Lepère et al. 2016, Fermani et al. 2021 ) and deep lakes with oxygenated hypolimnion (Mukherjee et al. 2017 ).All these diverse lacustrine ecosystems consistently reveal a substantial prevalence of unclassified sequences within numerous eukaryotic lineages.Compar ativ el y fe wer molecular biodiv ersity surv eys hav e been conducted in temperate lake environments (Lefranc et al. 2005, Boenigk et al. 2018, Mitsi et al. 2023 ).The water mixing patterns in holomictic freshwater lakes, where the water column is mixed in some seasons and remain stratified in other seasons, results in the recurring microbial communities' assembly processes.Deep dimictic lakes undergo mixing only during the spring and autumn months, maintaining stratification throughout the summer, and winter (Kirillin and Shatw ell 2016 ).Ho w e v er, climate c hanges influence lakes' mixing regimes (Adrian et al. 2009 ), which might pr ofoundl y impact these ecosystems by either enhancing or impeding vertical nutrient and dissolved gas fluxes (Råman Vinnå et al. 2021 ).Consequentl y, temper ate lakes, c har acterized by their water mixing patterns, offer a valuable opportunity to study seasonal pr otists' comm unity d ynamics (Le père et al. 2010, Medinger et al. 2010, Nolte et al. 2010, Mukherjee et al. 2017 ).
The Plankton Ecology Group model (Sommer et al. 1986(Sommer et al. , 2012 ) provides the best framework for describing the seasonal succession of phytoplankton and zooplankton in aquatic ecosystems.Se v er al studies (K ent et al. 2007, P av er et al. 2015, Woodhouse et al. 2016, Bock et al. 2018 ) have shown consistent temporal dynamics between eukaryotic phytoplankton and bacteria.Howe v er, prokaryotic and protist communities may show different temporal patterns over the course of the season.These differences could be due to variations in small-scale temporal patterns where prokaryotes and protists synchronize, as opposed to large-scale patterns where synchrony decreases due to changes in environmental conditions (Tammert et al. 2015, Obertegger et al. 2019 ).
Decades of r esearc h hav e unv eiled the pivotal role of physical factors, such as light, temper atur e, and turbulence in shaping of the microbial communities (Margalef 1978, Barton et al. 2015 ).Ho w e v er, seasonal succession is also governed by biological factors including organismal interactions (Drake 1990, Dakos et al. 2009, Logares et al. 2018, Bock et al. 2020 ).The advent of high-throughput DNA sequencing technologies has significantly bolstered our capacity to delineate microbial diversity and discern seasonal fluctuations within aquatic environments (Bunse and Pinhassi 2017, Giner et al. 2019, Grossart et al. 2020 ).Identifying these temporal patterns and determining their principal envir onmental driv ers ar e essential to r e v ealing the mec hanisms governing species succession and shaping community composition.Mor eov er, suc h inv estigations pr ovide v aluable insights into how climate change might alter these processes (Edw ar ds and Richar dson 2004, Siano et al. 2021, Caracciolo et al. 2022 ).
In this study, we conducted a metabarcoding investigation of the prokaryotic and protist communities in the temperate dimictic lake .We in v estigated small pr otists (size fr action 3-12 μm) and free-living bacteria (size fraction 0.2-3 μm) during the ice-free season to determine the dynamics of the community composition under pronounced seasonal gradients and to identify the main drivers of the communities in different seasons .We also in vestigated the influence exerted by abiotic parameters , e .g. temperature , organic carbon, and nutrient avail-ability, as well as biotic parameters on the microbial community structure.

Site description
Lake Ro ś (area: 18.08 km 2 ; maximum depth: 31.8 m, and mean depth: 8.1 m) is a meso/eutrophic glacial lake situated in northeastern Poland in the area of Masurian Lake District (53 • 38 -53 • 41 N 21 • 48 -21 • 59 E).It is a temperate dimictic lake, with biannual (spring and autumn) mixing e v ents.During the summer, Lake Ro ś experiences thermal stratification, leading to a pr onounced v ertical gr adient of dissolv ed oxygen (DO), r anging from oxic conditions in the epilimnion to near anoxia in the hypolimnion (Dawidowicz 1990 ).The lake periodically freezes during the winter months.Lake Ro ś has been a focal point for extensive ecological investigations throughout the 20th century, with particular attention given to macroph ytes, ph ytoplankton, zooplankton, macr oinv ertebr ates, and fish (e.g.Dawidowicz 1990, Jasser 1995, Pieczy ńska et al. 1998 ).The lake consists of two basins connected by a r elativ el y narr o w and shallo w channel (Fig. 1 A).The main southern basin is deeper (with a maximum depth of 31.8 m), maintains thermal stratification throughout the summer, and experiences common oxygen deficits within the hypolimnion.The second, northern basin is shallo w er (with a maximum depth of 9.3 m), and is pr edominantl y cov er ed by submer ged macr ophytes.This basin fr equentl y experiences summer destratification events, leading to complete mixing of its waters.

Sample collection
The sampling was conducted eight times during the ice-free season (from April to November) of 2019 in three sites within the lake (site A, B, and C on Fig. 1 A), which differ in their maximum depth r anging fr om 8 to 27 m (Fig. 1 B; Supplementary Table S1 ).The site A is located in the sheltered bay, in the vicinity of the deepest spot of the northern basin of the lake.Site B is located in the main basin of the lake, in the vicinity of the deepest spot of the lake.Site C is located on the periphery of the main basin, near the inflo w of tw o rivers: Świ ęcek and Konopka.In each site we obtained two samples of a total volume of 2 l using modified Bernatowicz sampler: one from the surface layer representing euphotic zone (3 m across all sites-A1, B1, and C1), while the second sample was taken from a depth of 2 m above the lake bottom (A2-6 m, B2-25 m, and C2-10 m).Samples were immediately filtered with a 150-μm plankton net to r emov e lar ge particles and m ulticellular or ganisms.Further filtr ation has been done sequentially with minimal (up to 200 hPa above atmospheric pressure) air pressure using Nucleopore filters (Whatman, Maidstone, UK) with 12, 3, and 0.2 μm pore size .T his process continued until filter clogging was detected, allowing us to obtain size fractions of 3-12 μm (r epr esenting small protists) and 0.2-3 μm (r epr esenting fr ee-living bacteria).Filters were then securely stored in −80 • C until the DNA extraction was conducted.
Planktonic animals filtered out from the samples using 150 μm plankton net were immediately preserved with 4% formalin.Subsequently, these specimens underwent thorough examination using dissecting microscopy.A subsample (10%-100% of the total sample volume, depending on zooplankton abundance) was taken and analysed to assess final abundance and taxonomic composition.Each planktonic animal in a subsample was identified and counted.Cladocer ans wer e identified to the species le v el, while copepods were identified at the order le v el, with nauplius (larv a) stages distinguished as a separate category (Rybak and Błędzki 2010 ).The densities of zooplankton taxa expressed in ind L-1 were calculated as the ratio of the number of animals observed within the subsamples to the corresponding total sample volumes.
Meteorological data for the lake area such as rainfall, and air temper atur e was deriv ed fr om publicl y av ailable r esources of the Institute of Meteorology and Water Management, Poland.Water temper atur e and oxygen le v el wer e measur ed during sampling at each sampling spot (A, B, and C) using a YSI Pr oODO m ultipar ametric probe, while the depth of the photic zone was measured using a portable light metre (LiCor LI-250A with spherical underwater PAR quantum sensor LI-193R) and the Secchi disc.The water was collected at two depths at each sampling spots and analysed for concentrations of biogenic (carbon, nitrogen, and phosphorus) and other (iron, manganese, and silicate) compounds by an external company (Wessling SA, Poland) ( Supplementary Table S2 ).

DN A extr action, DN A amplifica tions, and sequencing
For each sampling event and filter size (0.2 μm for prokaryotes and 3 μm for eukaryotic fr action), DNA was extr acted fr om one-fourth or half of the filter using the GeneMATRIX Soil DNA Purification Kit (EURx) for a total of 96 samples, according to manufacturer protocol including a single freezing/thawing step at −80 • C. Extracted environmental DN A w as quantified using NanoDrop (Thermo Scientific) and diluted to a concentration of 5 ng/ μl.Prokaryotic V4 hyperv ariable r egion of 16S rRN A gene (rDN A) w as amplified with Phusion High-Fidelity DNA pol ymer ase (ThermoFisher) using univ ersal pr okaryotic primers 515F-806R (Ca por aso et al. 2011 ) with further modifications (P ar ada et al. 2016 ) using r ecommended thermocycler conditions with 35 cycles .T he universal eukaryotic barcode V9 region of 18S rRNA gene (rDNA) was amplified using 1389F and 1510R primers (Amaral-Zettler et al. 2009 ), under recommended thermocycler conditions with reduced number of cycles (25) (de Vargas et al. 2015 ).All amplifications were done in triplicate in order to balance the variance within samples while obtaining adequate amounts of amplified DNA, combined, and then purified using a PCR clean-up kit (Syngen).The final concentration and quality of amplicons wer e a gain assessed by Nan-oDrop, and the library preparation and sequencing experiment on the Illumina MiSeq platform was performed by an external company (Genomed SA).The sequencing yielded 300 pair ed-end r eads targeted for 100 000 reads per amplicon.

Sequence analysis
Sequence quality c hec ks wer e conducted on r aw sequence data using FastQC (Andr e ws 2010 ), then sequencing ada pters wer e trimmed by trimmomatic (ILLUMINACLIP function) (Bolger et al. 2014 ).Subsequentl y, pr ocessed r eads wer e imported into the qi-ime2 environment and sequencing primers were removed using cutadapt (Martin 2011, Bolyen et al. 2019 ).Finally, D AD A2 denoising (minim um ov erla p = 12; max number of mismatc hes = 0 and consensus as a method for c himer a r emov al) was done for eac h sequencing batc h independentl y in qiime2 using dada2 denoisepaired function, after which amplicon sequencing variant (ASV) tables and feature tables were merged (Callahan et al. 2016 ).Taxonomic assignment of V4 16S rDNA was done using an RDP classifier encapsulated in qiime2 against SILVA99 138 database (Wang et al. 2007, Quast et al. 2012 ), and ASVs classified as 'Eukaryota', 'c hlor oplast', or 'mitoc hondria' wer e discarded.The assignment of V9 18S rDNA ASVs w as done b y USEARCH global alignment implemented in VSEARCH (Rognes et al. 2016 ) (minimum identity 60% and minimum query coverage 90%) against Protist Ribosomal Database PR2 4.14 (Guillou et al. 2012 ), pr epar ed following Tara Oceans guidelines (de Vargas et al. 2015 ).ASVs with the closest hit to a eukaryote, but with an identity lower than 80%, were assigned as an 'unknown eukaryote', and the r est wer e assigned to the best hit.In addition, prokaryotic V9 sequences were classified with usearch global alignment against the SILVA99 138 database (Quast et al. 2012 ).Befor e the main anal ysis, ASVs annotated as Metazoa, Embryophyta, Bacteria, or Archaea were also filtered out.Furthermor e, we assigned pr otistan ASVs to one of thr ee tr ophic gr oups ('phototr ophic', 'consumer', and 'par asitic') based on their taxonomic assignment and the published guidelines (Singer et al. 2021 ).Additionally, selected ASVs were manually annotated by liter atur e r esearc h and assigned to the gr oup 'mixotr ophic'.

Sta tistical anal yses and da ta visualiza tion
Statistical anal yses, suc h as alpha-div ersity, beta-div ersity, metaNMSD, ADONIS, ANOVA, Variation Partitioning Analysis or envifit and ancillary data visualizations, were performed in the R envir onment (v ersion 4.0.3)within RStudio IDE (Allaire 2012 ) using pac ka ges: v egan (Oksanen et al. 2023 ), qiime2R (Bisanz 2018 ), ggplot2 (Wickham 2016 ), phyloseq (McMurdie andHolmes 2013 ), a pe (P ar adis and Schliep 2019 ), and microbiome (Lahti and Shetty 2017 ).To r emov e the effect of inequality of sequencing depths, datasets were normalized using scaling with ranked subsampling-SRS (Beule and Karlovsky 2020 ).Eukaryotic and pr okaryotic pr e v alence anal yses wer e performed using the micr obiome pac ka ge (Lahti and Shetty 2017 ).The anal ysis was run for sampling sites A1, A2, B1, C1, and C2 (82 samples in total; sample B2 was analysed separately) from the whole sampling season.This division of the samples for further anal ysis r esulted fr om the fact that the hypolimnion was only formed in the sampling point B2, while all other points from deeper sampling points (A2 and C2) do not r epr esent the hypolimnion according to this definition, as no stable thermocline was formed ( Supplementary Fig. S1 ).
To compare these datasets, equal ranges of abundance and pr e v alence thr esholds wer e set and visualized using ggplot2.To inv estigate sync hr on y between 18S V9 and 16S V4 rRNA datasets, we used pca-based coinertia analysis implemented in the ade4 pac ka ge (Chessel et al. 2004 ).Two separate analyses were performed: (i) for samples r epr esenting epilimnion (sites A1, A2, B1, C1, and C2), and (ii) for the sample from hypolimnion (site B2).Each data table was first SRS normalized to get an e v en depth for each sample and Hellinger transformed according to Obertegger et al. ( 2019 ).The statistical significance of those analyses was c hec ked with the Monte-Carlo method implemented in RV.test from ade4 with 99 permutations .T he PC A-based coinertia analysis was visualized using ggplot2.To distinguish dead cells coming from upper layers of the lake fr om potentiall y living and thriving pr otistan linea ges in anoxic conditions, we compar ed r elativ e abundances of ASVs in epilimnion and hypolimnion layers during the stratification period (June-August), and only ASVs that ac hie v ed mor e than 2% of r elativ e abundance in at least one time point over this period.

Feature selection by Random Forests and Analysis of Compositions of Microbiomes with Bias Correction
To identify 16S rDNA and 18S rDNA ASVs whose abundance corresponded with seasons (two samplings in the spring, n = 10; three samplings in the summer, n = 15; and three samplings in the autumn, n = 15 per sequencing marker) in samples A1, A2, B1, C1, and C2 (epilimnion and metalimnion), we used supervised machine learning algorithm-Random Forests (RF) (Breiman 2001 ).This type of method has been pr ov en to perform well in classification of amplicon data (Hermans et al. 2020, Fang et al. 2022 ).To account for the substantial environmental variability in deeper layers associated with lake mixing, we excluded samples from the deepest point (B2), then ASVs that hav e mor e than 0.1% contribution were k e pt.ASV tables wer e then r enormalized after the filtration step and transformed using the scale function into scoring units .T he data was used to train RF models and then, Out-Of-Bag error was estimated.For each dataset, we picked 30 ASVs with the highest mean decrease Gini coefficient index scores, which corresponded to the highest impact on the classification of the samples, and visualized them with heatmaps and ordination plots.To confirm r esults fr om RF, w e emplo y ed the Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC) (Lin and Peddada 2020 ).The ANCOM-BC was run on absolute counts of the same samples as RF analysis with the option 'conserve', as recommended for the low number of samples, and P -values were adjusted using the Bonferroni correction.Subsequently, ASVs which wer e significantl y differ ent ( P -v alue < .05) in abundance wer e compared to ASVs selected by RF.For the visualization, we also added the 30 most abundant ASVs (estimated based on the sum of reads) for each of the analysed datasets.

Protist and bacteria di v ersity in Lake Ro ś
To investigate the plankton diversity in Lake Ro ś, we employed V9 18S rDNA and V4 16S rDNA amplicon sequencing of microbial community samples .T he samples were collected eight times over a 7-month period, ranging from 13th April to 18th November to r epr esent all the changes occurring during the vegetation season.We collected two size fractions of micr oor ganisms: the pr okaryotic fr action (0.2-3 μm) and small pr otists fr action (3-12 μm), fr om thr ee differ ent sampling locations and two depths, resulting in 96 samples (48 samples per molecular marker).For V9 18S rDN A, w e generated a total of 7628 535 reads (between ∼95 000 and 300 000 per sample) and inferred 7296 ASVs.Similarly, for V4 16S rDN A, w e obtained 5456 941 reads (between ∼40 000 and 170 000 per sample) and inferred in total of 6096 ASVs ( Supplementary Table S3 ).The r ar efaction curv e visualizations for both datasets confirmed that all samples r eac hed the plateau phase ( Supplementary Fig. S2 ).For eukaryotic amplicon analysis, we filtered out ASVs assigned to prokaryotes and Metazoa, resulting in a final dataset comprising 5191 ASVs accounting for 71% of the initial ASV count, with each sample containing ∼200-900 ASVs.Up to 1% of the prokaryotic sequences filtered out from the V9 18S rDNA datasets w ere classified as Ar chaea, including linea ges suc h as Methanosar cina, Nanosalinia, and Nanoar chaeia ( Supplementary Fig. S3A ).Two super gr oups ( sensu Burki et al. 2019 ) were highly prevalent across all samples: Cryptista (mainly cryptophytes-24.6%and katablepharids-3.4%)and TSAR super gr oup (telonemids-0.6%,str amenopiles-21.4%,alv eolates-24.5%, and rhizarians-6.7%)(Fig. 1 C).Regarding the V4 16S rDN A dataset, w e excluded eukaryotic sequences (pr edominantl y c hlor oplastic and mitochondrial) resulting in a dataset comprising 5690 ASVs, r epr esenting 93% of the initial count.Each sample contained around 300-800 ASVs , and o ver 99% of these were classified as Bacteria ( Supplementary Fig. S3B ).At the phylum le v el, the dominant bacterial phyla were Proteobacteria (40.8%),Bacteroidetes (25.2%), and Actinobacteria (24.5%) (Fig. 1 D).We also observed the temporal changes in the taxonomic composition of the protist ( Supplementary Fig. S4A ) and bacterial ( Supplementary Fig. S4B ) communities at different sites and depths.
To identify the cosmopolitan protists and bacteria in the surface layer (epilimnion) we applied the prevalence analysis.Only se v er al eukaryotic ASVs (16) were highly prevalent (occurred in more than 70% samples), and ther efor e could be considered as 'cor e' micr obiome for the ice-fr ee season ( Supplementary Fig. S5E ).Those belonged to Hacrobia (7), Alveolata (5), and Stramenopila (4), and accounted, on av er a ge, for 31% of r elativ e abundance of all pr otists.Onl y two ASVs were present in every analysed sample-ASV1 (classified to the genus Cryptomonas ) and ASV5 (Katablepharidophyta) ( Supplementary Fig. S5G ).In contr ast, the cor e pr okaryotic micr obiome was m uc h mor e numerous (50 ASVs) and accounted, on av er a ge, for 54% of the r elativ e abundance of prokaryotes ( Supplementary Fig. S5D and H ).
When analysing the core communities of various sites throughout the sampling season, we found that ∼40% of the core prokaryotic and eukaryotic ASVs wer e pr esent at each sampling site.Furthermore, 41.8% of eukaryotic ASVs and 34.4% of prokaryotic ASVs were only found at a single site .T he highest number of endemic eukaryotic ASVs were noted in site A (30.7%) with only 14.5% endemic prokaryotic ASVs.On other hand, the highest percentage of prokaryotic endemic ASVs was noted in the site C (15.6%) with only 3% of unique eukaryotic ASVs ( Supplementary Fig. S5A and B ). Mor eov er, the r atio between the mean and the maximum relative abundance of protist ('division' le v el) and prokaryotic ('class' level) taxa were much higher for eukaryotes (maximum 19-fold, noted for Discoba) than prokaryotes (maximum 7-fold) ( Supplementary Tables S4 and S5 ).Such differences suggest higher variability of the abundance of protists than pr okaryotes ov er the v egetation season.That could be also corr obor ated by the analysis of 'core' ASVs linked to each season ( Supplementary Fig. S5C and D ).Only 14.5% of eukaryotic core ASVs were consistently present in all seasons, comparing to 25% of prokaryotic ASVs .T he most unique core ASVs were identified in the summer for eukaryotes ( ∼40%) and autumn for prokaryotes (19%).

Protist and bacterial community structure is shaped by seasonal changes
We observed significant fluctuations in environmental parameters, including water temper atur e, oxygen le v els, light penetration, and the availability of chemical compounds across our sampling e v ents .T he surface water temper atur e v aried between 6 • C and 9 • C in April and November, while in June and August, it peaked at 23 • C. Notabl y, the r ecorded v ertical pr ofile of temper atur e and oxygen concentration during the period from June to August clearly indicated the occurrence of stratification ( Supplementary Fig. S1 ) in site B2.During the summer stratification phase, the deepest sampling location (B2-25 m) was below the pronounced oxycline and thermocline, and was characterized by the low temper atur e ( ∼10 • C), oxygen deficits and high concentrations of both organic and inorganic compounds ( Supplementary Table S2 ).Ho w e v er, in site A2 and C2 we only observed depletion of the oxygen in the summer months without clear thermocline, with the exception of Jul y, wher e in both sites we detected a homogeneous le v el of DO at all studied depths.
Observ ed c hanges of the envir onmental par ameters hav e discernible effects on the composition of microbial communities.In order to assess the diversity and similarity of these comm unities acr oss v arious sampling sites and time points, we conducted an analysis of α (richness) and β diversity.Shannon metrics varied between ∼2.3 and ∼5.4 for eukaryotes and between ∼3.5 and 5 for pr okaryotes.Notabl y, if we consider the sampling points, the highest eukaryotic diversity was observed in July at site B, while the lo w est occurred in May at the same site (Fig. 1 E).For the pr okaryotic comm unity, peak of div ersity was noted in August at site C, with the lo w est div ersity observ ed in September at site B (Fig. 1 F).Alpha diversity for eukaryotic plankton exhibited temporal dynamics, with lo w er diversity during spring and late autumn, ranging between ∼3 and 4, and higher diversity during summer, ranging between 4 and 5.This trend was supported by Anov a anal ysis, r e v ealing statisticall y significant differences between the 'month' and 'season' categories ( P -value < .001).Mor eov er, the differ ences in eukaryotic α div ersity wer e significant in the geogr a phical micr oscale and between depths within each month ( Supplementary Data 1 ).The pattern of α div ersity for pr okary otes w as quite different, ho w ever, it w as also fluctuating-after high diversity in April ( ∼4.5), it dropped during Ma y, J une , and J ul y (v aried between 4 and 4.2) and incr eased again in late summer and autumn ( ∼4.7).Anov a anal ysis indicated statistically significant differences between months, seasons , sampling sites , and depths within each month ( P -value < .05)for prokaryotic fraction ( Supplementary Data 2 ).To further understand the diversity of analysed prokaryotic and protist communities , we in vestigated their β diversity using the Unweighted UniFrac distance metric in conjunction with multidimensional scaling (MDS).In both the prokaryotic and protist fractions, we observed that sampling points formed three distinct clusters corresponding to the seasons (spring, summer, and autumn) points (Fig. 1 G and H), whic h was further confirmed with adonis anal ysis ( P -value < .001)and beta-dispersion analysis ( P -value > .05).Note worthy, in our anal ysis of β div ersity, we did not detect an y statisticall y significant differ ences between sampling sites for either 18S rDNA or 16S rDNA datasets ( Supplementary Data 1 and 2 ).Due to the formation of hypolimnion at site B2 from June to August ( Supplementary Fig. S1B ), the samples from this spot were consider ed separ atel y for further anal ysis.
Despite differences in diversity metrics and the sizes of core communities, the Principal Component Analysis-based Canonical Integr ation Anal ysis (PCA-based CIA) indicated that pr okaryotic and eukaryotic datasets displayed a coherent pattern of changes across seasons .T he analyses ha v e shown sync hr on y (RV = 0.87, P -value < .05)for shallo w er lay ers r epr esenting epi-and metalimnion (A1, A2, B1, C1, and C2) and for the sample B2 r epr esenting hypolimnion (RV = 0.96, P -value < .05)Ho w ever, the synchrony in the upper layers was disturbed in se v er al cases, especiall y during the summer, when it had a strong variation between prokaryotic and eukaryotic datasets compared to other seasons, which could cause the decrease of RV score ( Supplementary Data 3 , Supplementary Fig. S6 ).

RF analysis unveiled pivotal ASVs for seasonal community structures
To uncover seasonal changes in the eukaryotic and prokaryotic communities, we first focused on the 30 most dominant ASVs fr om eac h fr action (Fig. 2 ).Ho w e v er, of the dominant eukaryotic ASVs, most were persistent throughout the year, such as mixotrophic cryptophytes (ASV1, ASV6, ASV20, and ASV26) and chrysophytes (ASV10), as well as to the predatory katablepharidophytes (ASV5 and ASV22), ciliates (ASV24 and ASV28), and cercozoans (ASV27).Similar to eukaryotic plankton, most bac-terial ASVs persisted throughout the sampling period, but we could identify differences in their abundance between seasons.The dominant ASVs belonged to three phyla-Actinobacteriota (11), Bacteroidota (9), and Proteobacteria (10).
Since the dominant ASVs analysis did not explain well observed dynamic changes in microbial community structures across different seasons, we employed supervised machine learning (RF) and statistical (ANCOM-BC) methods, to identify ASVs significantly contributing to the shifts in community structure between seasons.For each dataset, we focused on the top 30 ASVs deemed most significant by the RF (RF-selected) model, as determined by the mean decrease Gini coefficient indexes ( Supplementary Fig. S7 ).Through the implementation of RF models, we effectiv el y categorized our samples into three distinct seasons [with an out-ofbag (OOB) estimate error rate of 0% for the 18S rDNA dataset and 2.5% for the 16S rDNA dataset].The ASVs r e v ealed selected by RF were also identified as statistically significant by the ANCOM-BC models, ho w e v er, onl y fiv e ASVs for eukary otes and tw o ASVs for prokaryotes exhibited overlap with the 30 most abundant ASVs ( Supplementary Fig. S8 ).
The 30 eukaryotic ASVs identified by RF analysis showed a nonuniform distribution of taxa across the sampling period, with ASVs clustered not only for the three seasons but also for the months (Fig. 3 A).Eac h of these clusters included ASVs r epr esenting differ ent tr ophic states, suc h as phototr ophs or mixotr ophs, consumers and par asites, whic h wer e assessed based on literatur e searc hes ( Supplementary Table S6 ).Spring was r epr esented by 12 RF-selected ASVs .T he presence of Bacillariophyta (ASV157) and parasitic fungi (ASV32)-chytrids-was consistent with the diatom bloom that typically occurs in April.During this period, we also observed a high abundance of photosynthetic Eustigmatophyceae (ASV46; Nannochloropsis ) and ASV59, which are assigned to predatory Stoeckeria (Dinoflagellata).In May, most ASVs belonged to the heter otr ophic assembla ge, whic h consisted of ASVs r epr esenting dinofla gellates (ASV9), ciliates (ASV31 and ASV85), cercozoans (ASV63, ASV68, and ASV113), and chrysophytes (ASV114), while parasites were represented by the genus Lagenidium (ASV164; Pseudofungi, Stramenopiles).In summer, the RF model and ANCOM-BC indicated the turnover of the observed taxa compared to the composition in spring.We were able to identify a more diversified group of photosynthetic or mixotrophic taxa belonging to c hlor ophytes (2), woloszynskioid dinoflagellates (2), Raphidophyceae (1), and Cryptophyta (1).Heter otr ophs were also diverse, including Ciliophora (3), Stramenopila-MAST-12 (1), Centroheliozoa (1), and Choanoflagellata (1).An ASV repr esenting par asites was also observ ed (ASV16; Perkinsea).Mor eover, taxa associated with summer were not evenly distributed, with a clear shift in the main ASVs assigned as primary producers and consumers between months .Among phototrophs , c hlor ophytes (ASV45; ASV100) r eac hed high r elativ e abundance in J une , woloszynskioid dinoflagellates (ASV133; ASV187) together with ASV111 (Raphidophyceae) in July, and ASV4 (Cryptophyta) in August.A similar succession was observed in the consumers: ciliate (ASV49), c hoanofla gellate (ASV80), MAST-12 (ASV65), and Chr ompodellida (ASV149) wer e particularl y conspicuous in J une , a centroheliozoan (ASV130) and a ciliate almost identical to Halteria (ASV124) in Jul y, while onl y a single ASVs (ASV85; Ciliophora) was noticeable in A ugust.A utumn can be roughly divided into two periods, the first (September and some spots in October) being c har acterized by a high r elativ e abundance of various Cercozoa (ASV12, ASV13, ASV68, and ASV113).In the second period (October and Nov ember), v arious consumers wer e pr esent, whic h wer e assigned to Cercozoa (ASV63), Str amenopila ( Labyrinthulea ; ASV132) and Amoebozoa (ASV126).Phototrophs, whic h wer e particularl y abundant in October, wer e assigned to the Bacillarioph yta (ASV33), Cryptoph yta (ASV4), and Eustigmatophyta (ASV46).Ho w e v er, some of the ASVs that wer e clearl y associated with a particular season were also present in other seasons .For example , the aforementioned ASV46, which was assigned to the Eustigmatophyta, was mainly present in early spring and late autumn, or the Cercozoa (ASV113), whose r elativ e abundance was high in either May or September.Grouping by season and individual months was also supported by the Bray-Curtis MDS visualization, whic h contr asts with the Br ay-Curtis MDS visualization of dominant ASVs, where data points were clustered together except for those from April ( Supplementary Fig. S9A  and C ). Analysis of the RF-selected ASVs within the prokaryotic dataset also r e v ealed a nonuniform distribution of taxa with four main groups of ASVs assigned to six phyla-Actinobacteriota (1), Bacteroidota (13), Cyanobacteria (1), Planctomycetota (1), Proteobacteria (11), and Verrumicrobiota (3) (Fig. 3 B).The first group of ASVs was associated with spring, with a considerable presence of ASVs assigned to the Bacteroidia (ASV132, ASV95, ASV142, and ASV154), two genera of Alphaproteobacteria-Sphingorhabdus (ASV36) and Rhodobacter (ASV117), ASV48 ( Pol ynucleobacter , Gamma pr oteobacteria), Verrucomicrobiae (ASV133), and CL500-3 clade of Phycisphaer ae (ASV55).Mor eov er, fiv e of them occurr ed mainl y in April (ASV55, ASV154, ASV36, ASV133, and ASV117), while four of them r eac hed higher abundance in May (ASV142, ASV132, ASV48, and ASV95).In summer, we observed two separate clusters of ASVsthe first in June and the second in July and August.The 'early summer' grouping consisted of six ASVs-classified as Bacteroidia, assigned to the env.OPS 17 group (ASV43 and ASV163), Verrucomicrobiae (ASV130 and ASV164), and Gamma pr oteobacteria (ASV37 and ASV179).The 'late summer' cluster of ASVs was formed by se v en ASVs belonging to Bacteroidia (ASV88, ASV131, ASV153, and ASV170), Cyanobium PC-6307 (Cyanobacteria, ASV139), and Gamma pr oteobacteria (ASV64 and ASV116).Inter estingl y, in a single sample in July (A1), we observed the recovery of the springassociated assemblage of ASVs.In autumn, ASVs were inconsistently distributed, with some ASVs being more abundant at the beginning (September and October) such as ASV159, assigned to Flavobacterium (Bacteroidia) and ASV57 ( Acinetobacter , Gamma pr oteobacteria) or at the end (October and November) such as members of Gamma pr oteobacteria (ASV56 and ASV57) and Bacteroidia (ASV63).Ho w ever, ASV45 ( Sphingobium , Alphaproteobacteria) and ASV96 ( Fluviicola , Bacteroidia) persisted throughout the whole autumn.In addition, there were ASVs that could not be explicitly associated with a specific time of the sampling period, such as ASV15 (Acidimicrobiia), which reached higher relative abundances in both spring and autumn, or ASV37, which was present in numerous samples in summer and autumn.The Bray-Curtis MDS visualization of the samples consisting of RF prokary otic ASVs sho w ed a slightl y differ ent arr angement to that of the eukaryotic ASVs, with three distinct clusters representing the seasons, albeit without a smooth transition between months.This is also in contrast to the Bray-Curtis MDS visualization of the pr okaryotic dominant ASVs, wher e all data points were clustered together ( Supplementary Fig. S9B and D ).

Distinct microbial community is established in the deep lake waters throughout the summer months
Throughout the summer months, spanning from June to August, Lake Ro ś undergoes a period of stratification, marked by the presence of a distinct thermocline that separates the epilimnion and hypolimnion la yers , eac h exhibiting markedl y differ ent envir onmental conditions.In comparison to the surface water, the hypolimnion layer was c har acterized with low temper atur e ( ∼10 • C) and the absence of sunlight and o xygen.Beta-di v ersity anal yses of prokaryotic and eukaryotic plankton r e v ealed that during this period, a discrete cluster of samples emerged within the hypolimnion layer at site B, located at a depth of 25 m (Fig. 4 A) implying the presence of a distinct microbial community.In contr ast, dissimilar micr obial comm unities at differ ent depths wer e not observed at other sampling sites, denoted as A and C (with deeper sampled fractions at 6 and 10 m), despite these locations also exhibiting episodes of anoxic conditions during this period ( Supplementary Table S2 ), suggesting that the distinct community of the site B is not only shaped by lack of oxygen.
An analysis of the taxonomic comparison of eukaryotic ASVs between the epilimnion and hypolimnion during the stratification period confirmed the differences in community structure ( Supplementary Fig. S3 ).Ho w e v er, r ather than observing a consistent community structure persisting throughout this period, we unexpectedly observed the formation of distinct assemblages for each of the three individual months (Fig. 4 ).In J une , five fungal ASVs ac hie v ed highest r elativ e abundances (45%), accompanied by a single diatom ASV assigned to Stephanodiscus (ASV3)-3% (Fig. 4 B; Supplementary Table S8 ).At the same time, only one ASV classified as Cryptophyta (ASV1) dominated the epilimnion (28%).In Jul y, onl y tw o ASVs w er e noted with high r elativ e abundancesthe pr e viousl y described diatom ASV3 (6%) and ASV261 (2.5%) annotated as c hoanofla gellate (Fig. 4 C; Supplementary Table S8 ).Fi-nally, the hypolimnion layer in August was mainly inhabited by a bodonid (ASV171)-9.5%, and Vermamoeba (ASV116)-3.5%,follo w ed b y an unkno wn eukary ote (ASV200)-7%, diatom (ASV3)-4%, and a r epr esentativ e of Katablepharidophyta (ASV5)-3.5% (Fig. 4 D; Supplementary Table S8 ).An analysis of the distribution of ASVs associated with the hypolimnion r e v ealed se v er al pr otistan ASVs ( Supplementary Table S8 ) across all samples, strongly implying their association with anoxic water conditions .T he exception was ASV116 (identified as Vermamoeba ), which was exclusiv el y found in site B2 at a depth of 25 m ( Supplementary Fig. S10A ).
The prokaryotic community structure was more uniform during the stratification period with a pr e v ailing pr esence of two phyla-Pr oteobacteria and Bacter oidota (Fig. 5 A-C).An anal ysis of the distribution of ASVs associated with the hypolimnion, at the taxonomic le v el of 'Famil y' acr oss all samples, r e v ealed a notabl y higher abundance of these ASVs in sites B2 and C2 in comparison to the other sampling locations ( Supplementary Fig. S10B ).

Abiotic and biotic factors influenced the microbial community composition
Through the incorporation of environmental parameters into our analyses, w e w ere able to discern the factors that exerted influence on the observed gradient or continuum of communities, as indicated by their correlation with the nonmetric multidimensional scaling (NMDS) axes (Fig. 6 A and B).In total, 18 environmental factors were tested using the envfit function, which entails fitting environmental vectors onto the NMDS ordination plot.This anal ysis r e v ealed that nine factors wer e significantl y corr elated with the ordination ( P -value < .01)for eukaryotic communities, and 15 factors exhibited significant correlations with bacterial communities .T he water and the air temper atur e, oxygen and Si concentr ation wer e the main factors shaping the structures of both communities.Other parameters such as light penetration, Secc hi disk visibility, concentr ations of NO 2 , NO 3 , NH 3 , Mn, TC, TN, F e , P, and the Trophic State Index were congruent with the pr okaryotic comm unity structur e. Dissolv ed or ganic matter (DOC) was found to be a significant driver exclusiv el y for the eukaryotic comm unity structur e.
The influence of abiotic factors on microbial community structur e was str onger; howe v er, at specific time points, biotic factors appear to assume a critical role .T his is in line with Variation Partitioning Analysis of eukaryotic and prokaryotic community structur es whic h sho w ed that 0.45 ± 0.02 of total variation w ere explained by abiotic factors, whereas 0.10 ± 0.02 of variation could be explained by biotic factors.Altogether these parameters explained more than 0.60 of total variation ( Supplementary Fig. S11 ).
Our focus centred on the potential grazing impact of zooplankton on the pr otist comm unity, and we conducted an analysis of the presence and abundance of the principal zooplankton groups.We used microscopic data collected for 11 zooplankton groups to perform abundance comparisons and correlation analyses and r e v eal their inter actions with pr okaryotic and eukaryotic plankton communities .T he abundance of zooplankton displayed temporal dynamics throughout the year, with the highest numbers occurring during the spring (May) and autumn (October) ( Supplementary Fig. S12 ).The marked increase in zooplankton abundance and pr esumabl y gr azing activity can explain an unexpected clustering pattern of protist groups across sampling sites within the same time points ( Supplementary Fig. S13 ).Such unexpected pattern was observed in the samples from May, where sample C1 exhibited a close relationship with B2, while sample C2 ASVs that r eac hed mor e than 2% of r elativ e abundance in a giv en sample wer e labelled.Eukaryotic gr oups ar e colour coded.clustered with sample B1.Within the C2-B1 cluster, a notably high r elativ e abundance of dinofla gellates (23%/31%) and ciliates (12%/30%) was observed, with a limited abundance of cryptophytes.In contrast, the B2-C1 cluster was dominated by cryptophytes (48%/70%), while ciliates (3%/13%) and dinoflagellates ( ∼2% in both cases) were less prevalent.The influence of zooplankton on planktonic protists was evident not only at the level of single ASVs but also at the community level.The analysis of the NMDS ordination plot with vectors fitting using envfit ( Pvalue < .01)implies that the larvae of Copepoda (nauplii), Bosmina longirostris , Chydorus sphaericus , Leptodora kindtii , and adult Cyclopoida.were the main drivers for shaping the community structure ( Supplementary Fig. S14 ).Spearman's correlation analysis sho w ed a more substantial impact of zooplankton on the protist comm unity compar ed to the pr okaryotic one, with 137 eukaryotic and 37 prokaryotic ASVs involved in statistically significant correlations ( P -value < .05 and r2 > | 0.4 | ) ( Supplementary Fig. S15 ).The analysis pointed to ciliates (33), c hlor ophytes (13), dinofla gellates (10), and cercozoans (10) as the taxa whose abundance correlates with zooplankton groups identified as Asplanchna sp. ( 44) and D. brachyurum (44).

Di v ersity of eukaryotic and prokaryotic plankton in Lake Ro ś
The 16S rDNA and 18S rDNA metabarcoding allowed us to identify the v ast div ersity of planktonic bacteria and small protist community of the Lake Ro ś.The prevalence of Cryptophyta and Katablepharidophyta (Hacr obia), Ciliophor a and Dinofla gellata (Alveolata), and Bacillariophyta (Stramenopila) among eukaryotes (Fig. 1 C), and Actinobacteria, Bacter oidota, and Pr oteobacteria among prokaryotes (Fig. 1 D), was in line with pr e vious r eports on the microbial community composition in temperate zone lakes (Liu et al. 2015, Cruaud et al. 2019, 2020, Kiersztyn et al. 2019 ).
We have also observed dynamic changes in the taxonomic composition of both prokaryotic and pr otistan comm unity acr oss the seasons, indicating a temporal succession of species .T his w ell-kno wn phenomenon finds support in both richness and betadiv ersity anal yses of the Lake Ro ś comm unity (Fig. 1 E-H).Ev en though the bacterial and eukaryotic communities in Lake Ro ś exhibited dynamic changes over time in different scale, the coinertia analysis based on PCA revealed their synchrony, which might be a result of direct biotic interactions or response to the same envir onmental factors (Boc k et al. 2018 ).In Lake Ro ś, sync hr on y was observed in all spots, ho w ever, in shallo w er w aters experienced disruptions on se v er al occasions, particularl y during the summer season.These disturbances could be due to factors such as the strong dominance of certain taxa, selective grazing of zooplankton and climatic disturbances such as heavy rainfall and additional mixing in the shallo w er parts of the lake ( Supplementary Fig. S6 , Supplementary Table S2 ).
While it is undebatable that the geogr a phical distance has impact on the pr otist comm unity structur e (Sc hiaffino et al. 2016, Boenigk et al. 2018 ), m uc h less attention was put into diversity within single water body.Most of the recent studies, whether focused on single or multiple lakes , ha v e c har acterized micr obial communities based on a single sampling site per lake (Simon et al. 2015b, Sieber et al. 2020 ), ho w e v er, the statisticall y significant effect of the sampling site was observed for larger fresh water bodies, such as Lake Baikal (David et al. 2021 ).Our sampling strategy encompassed three distinct sampling sites in Lake Ro ś, Figure 6.Impact of physicochemical factors on microbial community structure.Sample composition of (A) eukaryotic and (B) prokaryotic datasets based on NMDS analysis with envfit (vegan package) fitted physicochemical parameters ( P < .05).The colours refer to the seasons, the shapes to the sampling sites.
enabling the detection of variations in microbial community diversity within a single lake ( Supplementary Fig. 5 ).Although betadiv ersity anal ysis highlighted seasonality as the predominant factor shaping the communities, with the 'site' factor not attaining statistical significance, we observed disparities between the core taxa composition in surface layers of sites A, B, and C. Particularl y note w orthy w as the eukary otic cor e comm unity at site A, r epr esenting 30.7% of all ASVs .T he prokaryotic community exhibited lo w er variability, y et the most distinct communities w ere discerned at sites A and C, constituting 14.5% and 15.6%, respectiv el y.These findings suggest that the micr obial comm unity r esponse to local physicochemical and biological factors is significantly influenced by the hydrological characteristics of the lake.For instance, site A was the shallowest sampling site (6 m deep) experiencing consistent year-round mixing, separated from other points by a narrow channel.In contrast, site C is potentially impacted b y inflo ws of allochthonous matter carried by two small rivers that collect water from surrounding farmlands .T hese disparities at the microscale underscore the importance of carefully selecting sampling locations, as it can significantly impact the accurate identification of distinct microbial communities.
Micr obial comm unities c hange ov er m ultiple timescales (i.e. from hours , da ys , or weeks to seasons) in response to a multitude of abiotic and biotic factors (Fuhrman et al. 2015 ).In temperate lakes, the recurring seasonal patterns strongly influence both prokaryotic (Crump andHobbie 2005 , Rusak et al. 2009 ) and eukaryotic (Simon et al. 2015b ) community composition.Howe v er, these c hanges ar e often examined at the le v el of functional gr oups, encompassing primary pr oducers and consumers, or lar ge taxonomic groups.Our metabarcoding analysis, thanks to the application of machine learning and statistical methods , pro vides more detailed insights into the succession of protist and bacteria at the le v el of ASVs.

Protist communities of Lake Ro ś
The dominant eukaryotic ASVs during most of the ice-free season belonged mainly to the mixotrophic groups (especially cryptoph ytes and chrysoph ytes) and eukary ov orous katablepharids (Fig. 2 A; Supplementary Table S8 ).The continuous presence and predominance of these groups suggests that they do not have a significant influence on the formation of unique seasonal comm unities.Ne v ertheless, it is worth noting that the pr e viousl y mentioned dominant protistan plankton groups could potentially under go c hanges during winter, as pr e vious studies with annual sampling have shown (Cruaud et al. 2019 ).Mixotrophic plankton, particularl y pha go-mixotr ophic or ganisms that combine photosynthesis with pha gotr ophy, ar e of particular interest due to their dual role as both producers and consumers (Millette et al. 2023 ).On the one hand, their mixotrophic potential can be an advantage and increase their flexibility in adapting to changing environmental conditions by grazing on bacteria, and thus displacing phototrophic species (Selosse et al. 2017 ).On the other hand, increased nutrient and organic matter inputs to certain lakes may affect mixoplankton by promoting bacterial prey population growth and limiting light availability due to humic substance absorption in the water column (Wilken et al. 2018 ).
Emplo y ed RF models and Ancom-BC both highlighted the role of less abundant linea ges c har acteristic for eac h season (Fig. 3 ).In fact, the set of ASVs identified b y RF w as v ery differ ent fr om the dominant ones .T he RF model selected se v er al linea ges like phototr ophic and heter otr ophic dinofla gellates, phototr ophic eustigmatophytes and c hlor ophytes, or heter otr ophic cercozoans and ciliates, which might be important for the functioning of the micr obial comm unity.Importantl y, eac h season was r epr esented by different lineages of protists (different ASVs), even though the same functional categories such as phototrophs , consumers , and par asites wer e pr esent thr ough the whole sampling period (Fig. 3 A).Among primary producers, beside re presentati ves of expected gr oups suc h as diatoms and c hlor ophytes, we identified Eustigmatophyceae ( Nannoc hloropsis ), whic h wer e pr e viousl y documented in spring blooms of freshwater lak es (Fawle y and Fawley 2007 ).Through taxonomic classification and extensive literature sear ches, w e w ere able to classify some of the detected dinoflagellates as primary pr oducers, particularl y those associated with photosynthetic genera Asulcocephalium (ASV133) and Leiocephalium (ASV187) (Takahashi et al. 2015 ).The diversity of consumers, on the other hand, was significantly greater in all seasons.We were able to identify taxa that are possibly responsible for grazing within the eukaryotic community.This included two dinoflagellates related to Gyrodinium (ASV9) and Stoeckeria (ASV59), as well as a ciliate from the Balanion genus (ASV49) in the spring.Gyrodinium has been pr e viousl y r eported as a major grazer of diatoms in marine systems, as opposed to ciliates, which are less capable of consuming large prey (Saito et al. 2006 ).In May, three eukary ov orous cercozoans ( Supplementary Table S7 ) belonging to two closely r elated gener a-Protapis (ASV68) and Cryothecomonas (ASV13 and ASV63), known as typical marine diatom pr edators (Dr ebes et al. 1996, Schnepf and Kühn 2000, Moustaka-Gouni et al. 2016 ), may be involved in controlling the decline of diatom blooms.We propose that other RF-selected heter otr ophs wer e enga ged in gr azing on bacteria.In the spring, that might be ASV114 ( Pedospumella , Chrysophyceae) known to be important bacterioplankton predator in freshwater ecosystems (Šimek et al. 2013 ).Ho w e v er , during the summer, bacterivorous protists displayed greater taxonomic div ersity, including 'r ar e taxa' r epr esenting colpodellids, stramenopiles , heliozoans , and c hoanofla gellates .T his disco very of r ar e taxa underscor es their importance for the summer micr obial comm unity and str ongl y suggests their seasonality (Schiwitza et al. 2020, Za gumyonn yi et al. 2022 ).Ne v ertheless, further r esearc h is r equir ed to inv estigate their ecological r oles in lakes.Among the heter otr ophic or ganisms, we also identified potential decomposers, namely the ASV132, related to Thraustochytrium sp.(Labyrinthulomycetes, Stramenopila), which represents a significant group of marine and freshwater saprotrophic eukaryotes known for their ecological role as decomposers (Pan et al. 2017, Morabito et al. 2019, Xie et al. 2022 ).While we can only speculate on the exact role of this lineage in freshwater ecosystems, it likely contributes to the decomposition of biomass from ongoing summer c y anobacterial blooms .T he pr esence of c hytrids (ASV32), perkinsids (ASV16), and pseudofungi (ASV164 and ASV23), in the seasonal pr otist comm unities suggest their potential influence on the diversity and dynamics of freshwater ecosystems (Mangot et al. 2009 ), but also raises concerns for host-parasite interactions, which might be impacted with the increase of eutrophication (Budria 2017 ).Of particular note is the remarkable diversity of Perkinsea-related ASVs, comprising a total of 106 distinct ASVs .T his parasitic group, with a wide host range spanning from dinoflagellates to animals, poses a potential risk in freshwater envir onments, wher e it has been linked to the occurrence of mass mortalities among amphibian species (Isidoro-Ayza et al. 2017, Itoïz et al. 2022 ).The majority of RF-selected ASVs were exclusiv el y pr esent in a single season.Ho w e v er, ther e wer e certain taxa that r ecurr ed acr oss m ultiple seasons, impl ying a potential ada ptation to specific biotic and abiotic factors, such as the presence of pre y, n utrient availability, or favourable temperature conditions.For instance, the repeated presence of eukary ov orous dinoflagellates (ASV9; ASV59) could be linked to their specialization in preying upon diatoms, suggesting a specific ecological niche associated with diatom availability.

Bacterial communities in Lake Ro ś
Most of the dominant ASVS of the bacterioplankton community remained relatively consistent throughout the study period.The group of dominants comprised various re presentati ves of Actinobacteria, Acidimicr obiia, Gamma pr oteobacteria, and Bacteroidota (Fig. 2 B), previously described as a major component of freshwater ecosystems (Kiersztyn et al. 2019, Cruaud et al. 2020 ).The RF-based model and ANCOM-BC emphasized four major groups of prokaryotic ASVs associated with seasons, primar-ily affiliated with Proteobacteria, Bacteroidota, and Verrucomicrobiota (Fig. 3 B).These results suggest the division into dominants, whic h ar e widespr ead gener alists and less numer ous specialists.Inter estingl y, specialists in aquatic systems are frequently involved in the degradation of dissolved organic carbon (DOCp) fr om phytoplankton, whic h is pr oduced by exudation or cell lysis (Sarmento et al. 2016 ).Combined with the fact that different types of phytoplankton promote the growth of different bacterial groups (Sarmento and Gasol 2012 ), our results suggest that many of the observed intermittent occurrences of ASVs might be a result of such associations .T he ASVs selected b y RF w ere not onl y mor e div erse (e v en at the phylum le v el) than the dominants, but also belonged to groups previously reported to be involved in the degradation of certain organic compounds produced by phytoplankton.This is most evident in April, where along with diatom bloom we observed the ASV117 ( Rhodobacter ), ASV113 (Verrucomicrobiae), and members of the Bacteroidota (ASV132, ASV142, and ASV154), whic h ar e either specialized in the degradation of diatom DOCp or are generally associated with diatom blooms (Tada et al. 2012, Orellana et al. 2022 ).Furthermore, a similar assemblage of taxa (Verrucomicrobiae and Bacteroidota) was also present in J une , although we cannot explicitly point to the source of the or ganic matter, whic h could be either c hlor ophytes and r a phidophytes, or other algae (Fig. 3 A).In the 'late summer' assemblage, we also observed taxa such as the genus Fluviicola (ASV170), which has been reported to be closely associated with blooms of primary producers (Eckert et al. 2012 ) such as those of colonial species like Microcystis or Dolichospermum spp., which are typicall y observ ed during c y anobacterial blooms (Woodhouse et al. 2018 ).

Microbial community of the anoxic hypolimnion
The majority of r esearc h on hypolimnion ecology focused on deep freshwater lakes with oxygenated hypolimnion, leading to the identification of distinct micr obial comm unities in these environments (Okazaki and Nakano 2016, Mukherjee et al. 2017 ).Howe v er, anthr opogenic eutr ophication of lakes and climate change increases the number of lakes experiencing anoxic conditions in their hypolimnion during the summer months and it is reasonable to anticipate that such anoxic environments would also harbour unique micr obial comm unities.Sur prisingl y, studies inv estigating pr otistan comm unities in anoxic hypolimnion have been r elativ el y scarce thus far.A case in point is the Lake Ro ś, which exhibits anoxic conditions within its hypolimnion during the summer.T his en vir onmental c har acteristic is further r eflected in the establishment of a distinct hypolimnetic micr obial comm unity (Fig. 4 A).In suc h envir onments, the main driv ers of micr obial community is decomposition of organic matter by methanogens-a phenomenon well-documented in other lakes c har acterized by high redox potential during stratification (Reis et al. 2022, Shi et al. 2022, Steinsdóttir et al. 2022 ).While w e w ere not studying arc haeal comm unity, we did identify methanogenic archaea Methanosarcina ( Supplementary Fig. S3 ) amplified b y eukary otic V9 primers in the B2 samples (Choi andPark 2020 , Carr andBuan 2022 ).Consequently, we conclude that methane may be oxidized b y Meth ylobacter (ASV30) that uses alternate electron acceptors under anoxic conditions (Hao et al. 2022 ).Recent reports, also suggest the strong syntrophy between Methylobacter and Methylotenera (ASV7), which was the dominant methylotrophic genus that accounted for up to 9.5% of r elativ e abundance .T heir relationship couples nitrogen and carbon (C1) cycles with the extensive use of nitrates as alternative electron acceptors, which helps transfer carbon from methane to other members of hypolimnetic food web, such as Bacteroidota (ASV9 and ASV29), or the methylotrophic Methylophilus (ASV8) (Van Grinsven et al. 2021 ).In addition to the observed impact of biomass influx from the upper layer of the lake on the microbial community structure within the hypolimnion, there are reports suggesting the significance of other compounds often found in high concentr ations, suc h as iron and manganese, in shaping the community.High concentrations of ir on wer e observ ed in samples from site B of Lake Ro ś and can be associated with the growth of specific bacterial taxa, such as Candidatus Omnitr ophus whic h r elies on ir on for ma gnetosome biosynthesis (Kolinko et al. 2016 ).
Within this ecosystem, the majority of identified protists ar e bacterivor es .T hese include the genus Bodo (ASV171) and a c hoanofla gellate (ASV261), whic h hav e pr e viousl y been r eported in various marine and brackish anoxic environments (Bernard et al. 2000, Stock et al. 2009 ).Notably, we also identified an unknown eukaryotic lineage represented by ASV200 (with a low sequence identity ∼84% to an unknown eukaryote), suggesting the potential for the discovery of novel eukaryotic lineages within such ecosystems.Our data show that a well-established hypolimnetic community is formed in August, whereas, during the months of June and July, eukaryotes were dominated by ASVs related to Fungi or stramenopiles (diatoms), suggesting a potentially significant influence from the influx of dead algae, while eukaryotic ASVs associated with anoxia were present during this period in low abundance .T hese findings suggest that sinking dead cells and organic particles are vital early contributors to eukaryotic and bacterial plankton community development in the hypolimnion as a source of organic matter.An open question remains regarding the persistence of hypolimnion-associated lineages during the biannual water column mixing in spring and autumn.Most likely their refuge is in the sediments, often inhabited by obligately anoxic benthic microbiomes (Bernhard et al. 2014, Gomaa et al. 2022 ).

Abiotic and biotic factors shaping microbial community
We conducted an analysis of various abiotic factors pr e viousl y shown to influence micr obial comm unities (Boc k et al. 2020 ).Among these factors, only a subset was found to play role in shaping the microbial communities of Lake Ro ś.T emperature, DO , and silicon concentration emerged as the primary drivers of both comm unities' structur es (Fig. 6 ).While temper atur e and DO le vels have been recognized as significant factors shaping microbial community composition in many studies (Liu et al. 2013 ) Oliverio et al. 2018, Boenigk et al. 2018, Mikhailov et al. 2022, Shang et al. 2022 ), the silicon concentration seems to be mainl y r elated to temperate lakes, with spring and autumn mixing events (Panizzo et al. 2018, Kong et al. 2021 ).In Lake Ro ś, the silicon concentration r eac hed its maxim um during the autumn mixing, follo w ed b y the diatom bloom in spring, leading to a decrease in silicon concentration due to its utilization by diatoms for building their cell walls.A similar trend has been observed in Lake Baikal (Mikhailov et al. 2022 ).Other factors, including trophic status and phosphorus and nitr ogen concentr ations, corr esponded with the pr okaryotic comm unity structur e, but did not show str ong impact on the eukaryotic comm unity.These r esults align with pr e vious r esearc h in freshwater lakes, suggesting that changes in bacterial communities ar e mor e closel y linked to physicoc hemical patterns compared to protist communities (Bock et al. 2020 ).
Observed distinct seasonal patterns are not only more evident for eukaryotes than prokaryotes, but they are also more sta-ble when facing short-term environmental fluctuations (Jacobsen and Simonsen 1993, Stockwell et al. 2020 ).For example, in sample A1 fr om Jul y, we observ ed the occurr ence of ASVs associated with a spring bacterial assemblage .T his was pr obabl y due to a drop in water temper atur e fr om ar ound 23 • C in June to around 20 • C in July, caused by rainfall and low air temperatures .T his drop in temper atur e led to an increase in the r elativ e abundance of diatoms, particularly ASV3, by up to 5%, which most likely stimulated the growth of diatom DOCp specialists (Fig. 2 A and 3 B; Supplementary Table S2 ).This example, among others, clearly shows that the monthly sampling scheme is not sufficient to determine the universal factors influencing the changes in microbial communities, because their turnover occurs rather in days than weeks (Šimek et al. 2014 ).In our study, we also unveiled the significant influence of zooplankton on the diversity of protists ( Supplementary Figs S13 and S14 ), a critical component of freshw ater food w ebs, as top-do wn r egulators of lar ger pr otists (Lu and Weisse 2022 ).This impact became evident on a global scale thr ough the corr elation between pr otistan ASVs and zooplankton cell counts, particularly in the case of the predatory omnivorous r otifer Asplanc hna , and the cladocer an Diaphanosoma brac hyurum with pr efer ence for smaller ( < 3 μm) particles (Chang et al. 2010, Nandini et al. 2021 ).Furthermore, at a local scale, this influence was evident in beta-diversity patterns.We observed an unusual clustering pattern among samples from May ( Supplementary Fig. S13 ), suggesting that zooplankton, thr ough their gr azing activity (top-down selection) during clear water phase, could dr amaticall y alter the local composition of protists in a specific region within the lake.

Climate change impact on microbial communities in dimictic temper a te lakes
Our data point to the role of temperature and oxygen levels in the formation of different planktonic communities (Fig. 6 ).Rapid changes in these factors, driven by climate change, are expected to have increasingly profound effects on freshwater lake ecosystems, thereby impacting biodiversity and ecosystem functions.Notabl y, surface water temper atur es of fr eshwater ecosystems are rising at an accelerated rate compared to air and ocean temper atur es (O'Reill y et al. 2015, Dokulil et al. 2021 ).Recent studies hav e demonstr ated that warming pr edominantl y leads to a decr ease in fr eshwater plankton div ersity (Rasconi et al. 2017, Bergkemper et al. 2018, Verbeek et al. 2018, Da Silva et al. 2019, Celewicz and Gołdyn 2021 ).Ho w ever, the effects are multifaceted, often resulting in shifts in community structures, particularly tow ar ds green algae dominance (Rasconi et al. 2017, Yu et al. 2018, Zhang et al. 2021, Beng et al. 2023 ).Considering the complexity of these changes, it is essential to assess diversity at an a ppr opriate taxonomic le v el, since the shifts might not be evident at higher taxonomic le v els.Our study illustr ates that the abundance of specific ASVs can exhibit dynamic changes over time, even when the abundance of the broader taxonomic groups they belong to remains r elativ el y stable (Fig. 2 A and B).
Additionall y, climate c hange is driving a decline in DO le v els in aquatic ecosystems, affecting lakes, coastal zones, and open oceans globall y (Sc hmidtk o et al. 2017, Br eitbur g et al. 2018, Limburg et al. 2020, Jane et al. 2021 ).Lar ge-scale anal yses r e v eal that the majority of lakes (over 70%) are experiencing increases in o xygen-de pleted w ater (J ane et al. 2023 ).This trend is significant since reduced DO concentrations can be observed during late summer periods due to changes in stratification characteristics (Jane et al. 2023 ), including earlier onset of seasonal stratification and less frequent mixing events (Woolway and Merchant 2019 ).We have identified a distinct hypolimnetic community (Fig. 4 ; Supplementary Table S8 ) that thrives in oxygendepleted waters .T his community is primarily composed of kinetoplastids (ASV171), c hoanofla gellates (ASV261), an amoebozoan (ASV116), and ASV200 from a novel protist group.While their role in the lake's ecosystem is not as well-understood as that of epilimnetic plankton, the increase in hypoxic zones highlights the critical need to explore the taxonomic and functional diversity of this community.

Conclusions
In our metabarcoding study of a temperate dimictic Lake Ro ś, we gained insights into the taxonomic composition and community structure of small protist and free-living bacteria during the icefree period at the unprecedented level of single ASVs.Leveraging RF and ANCOM-BC analyses, we identified ecologically functional clusters of eukaryotic and prokaryotic ASVs that were associated with different seasons.Seasonal changes were mainly associated with consumer groups such as cercozoans, and parasitic taxa such as Pseudofungi and Chytridiomycota.In contrast, the generalist ASVs, at least during the ice-free season, were mainly phototr ophic and mixotr ophic or ganisms, suc h as Cryptophyta, and pr edators, suc h as Katablepharidophyta.The pr okaryotic diversity could also be divided into generalists (such as Actinobacteria) and specialists, which are a diversified group of taxa that are most likel y involv ed in r ecycling or ganic matter, suc h as DOCp, abundant at certain time points.Significant differ ences wer e also observ ed between micr obial comm unities in the epilimnion and hypolimnion, with k e y hypolimnion-specific taxa identified, including Choanoflagellata, Amoebozoa (Lobosa), Discoba (Kinetoplastida), and a putativ e nov el linea ge (ASV200).These taxa likel y feed on a prokaryotic community driven by organisms involved in C1 cycle, such as methanogens and methanotrophs .T he observ ed differ ential seasonal patterns in pr otistan and pr okaryotic communities align with their distinct responses to environmental factors.Eukaryotes exhibited differ ent r esponses compar ed to pr okaryotes, particularl y to the three main factors of temperature, oxygen, and silicon concentration.While these factors affected both groups , other en vironmental variables primarily influenced bacterial communities (bottom-up regulation).In contrast, zooplankton composition and abundance exerted a more pronounced top-down influence on the eukaryotic community compared to the prokaryotic community.

Ac kno wledgements
Sampling was carried out with the equipment of the Hydrobiological Field Station of the Faculty of Biology.We would like to thank e v eryone who helped us with the sampling, especially the head of the station Mirosław Ślusarczyk, and Maria Turłaj.

Figure 1 .
Figure 1.Sampling scheme and ov ervie w of the microbial community structure.Location of sampling sites A, B, and C in Lake Ro ś (A) with the sampling depths (m) for each sampling site (B).Tr eema ps r epr esent the ov er all composition of r elativ e abundances for 18S V9 rRNA at the 'Class' le v el (C) and 16S V4 rRNA at the 'Phylum' le v el (D).Boxplots illustrating Shannon alpha diversity for each month and depth (class 'Upper' represents samples A1, B1, and C1 and class 'Lo w er' r epr esents A2, B2, and C2) for 18S V9 rRN A (E) and 16S V4 rRN A (F) datasets, sites are colour coded.Ordination plots based on Unweighted UniFrac MDS for 18S V9 rRNA (G) and 16S V4 rRNA (H), with sampling sites marked with shapes and seasons marked with colours.

F igure 2 .
Dominant eukary otic and prokary otic ASVs in the epilimnion during the entir e sampling period.(A) Heatma p depicting the abundance of 30 dominant ASVs for eukaryotes and (B) heatmap depicting the abundance of 30 dominant ASVs for prokaryotes .T he trophic states for eukaryotes and phyla for prokaryotes are colour coded.

Figure 3 .
Figure 3. Model selected eukaryotic and prokaryotic ASVs in the epilimnion throughout the entire sampling period.(A) Heatmap depicting the abundance of 30 ASVs selected by the RF model for eukaryotes (B) Heatmap depicting the abundance of 30 ASVs selected by the RF model for prokaryotes .T he trophic states for eukaryotes and phyla for prokaryotes are colour coded.

Figure 4 .
Figure 4. Pr efer ences of eukaryotic taxa between the hypolimnion and the epilimnion during stratification at site B2 (25 m).(A) Dendrogram of the eukaryotic samples (18S V9 rRNA) based on the Unweighted UniFrac metric ('complete' clustering method), with the hypolimnion samples marked with red arrows.Comparisons of relative abundances (%) of ASVs between the surface (B1) and bottom part (B2) for June (B), July (C), and August (D).ASVs that r eac hed mor e than 2% of r elativ e abundance in a giv en sample wer e labelled.Eukaryotic gr oups ar e colour coded.

Figure 5 .
Figure 5. Pr efer ences of prokaryotic taxa between the hypolimnion and the epilimnion during stratification at site B2 (25 m).(A) Dendrogram of the prokaryotic samples (16S V4 rRNA) based on the Unweighted UniFrac metric ('complete' clustering method), with the hypolimnion samples marked with red arrows.Comparisons of relative abundances (%) of ASVs between the surface (B1) and bottom part (B2) for June (B), July (C), and August (D).ASVs that r eac hed mor e than 2% of r elativ e abundance in a giv en sample wer e labelled.Pr okaryotic classes ar e colour coded.